A short description of the post.
Important to set fig.retina to be 3 in order for subsequent visualization to be sharp.
packages <- c('maptools','sf','raster','spatstat','tmap','tidyverse', 'plotly','ggthemes')
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
}
library(p, character.only=T)
The code chunk import shapefile using st_read() of sf package. The output object is in tibble sf object class.
Pgl_sf and mpsz_sf geospatial data sets into RStudio.
Note: - Do not paste directory in dsn
Pgl_sf <- st_read(dsn="data/shapefile", layer="Punggol_St")
Reading layer `Punggol_St' from data source
`C:\Users\User\Desktop\IS415\lye-jia-wei\IS415_blog\_posts\2021-09-13-in-class-exercise-5\data\shapefile'
using driver `ESRI Shapefile'
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM
mpsz_sf <- st_read(dsn="data/shapefile",layer="MP14_SUBZONE_WEB_PL")
Reading layer `MP14_SUBZONE_WEB_PL' from data source
`C:\Users\User\Desktop\IS415\lye-jia-wei\IS415_blog\_posts\2021-09-13-in-class-exercise-5\data\shapefile'
using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
Projection is SVY21
read_rds() of readr package is used instead of readRDS() of base R is used. This is because output of read_rds() is in tibble object.
Note:
No need to set path
Note that there are some data issue in childcare data frame with the Lat and Lng column, the columns should be in numeric data type. The coordinate fields are in decimal degrees. Hence, wgs referencing system is assumed.
childcare <- read_rds("data/rds/childcare.rds")
CHAS <- read_rds("data/rds/CHAS.rds")
Note:
st_as_sf() accept coordinates in character data type.
No significant difference whether coordinate data type is change
Alpha is to set transparency
tmap_mode("view")
tm_shape(childcare_sf)+tm_dots(alpha= 0.4, col="blue", size=0.05)+ tm_shape(CHAS_sf)+tm_dots(alpha= 0.4, col="red", size=0.05)
Note:
as_Spatial() of sf package.
childcare <- as_Spatial(childcare_sf)
CHAS <- as_Spatial(CHAS_sf)
mpsz <- as_Spatial(mpsz_sf)
as.SpatialPoint() of as.SpatialPolygon() of maptools package
childcare_sp <- as(childcare, "SpatialPoints")
CHAS_sp <- as(CHAS, "SpatialPoints")
MPSZ_sp <- as(mpsz, "SpatialPolygons")
This is to drop the geographical projection information using as.ppp() of maptools package.
childcare_ppp <- as(childcare_sp, "ppp")
CHAS_ppp <- as(CHAS_sp, "ppp")
childcare_ppp_jit <- rjitter(childcare_ppp, retry = TRUE, nsim=1, drop = TRUE)
any(duplicated(childcare_ppp_jit))
[1] FALSE
Output return a ‘False’ which suggest that there are no more duplicate.
Note:
CHAS_ppp_jit <- rjitter(CHAS_ppp, retry = TRUE, nsim=1, drop = TRUE)
any(duplicated(CHAS_ppp_jit))
[1] FALSE
Output return a ‘False’ which suggest that there are no more duplicate.
Note:
Extracting planning sub zone PUNGGOL from data and then planning area.
pg <- mpsz[mpsz@data$PLN_AREA_N=="PUNGGOL",]
Note:
PUNGGOL planning areapg_sp <- as(pg, "SpatialPolygons")
Convert into owin object which is a list of length 5.
pg_owin <- as(pg_sp, "owin")
childcare_pg <- childcare_ppp_jit[pg_owin]
CHAS_pg <- CHAS_ppp_jit[pg_owin]
plot(childcare_pg)

L_childcare <- envelope(childcare_pg,
Lest,
nsim = 99,
rank = 1,
gloval = TRUE)
Generating 99 simulations of CSR ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
Done.
title <- "Pairwise Distance: L function"
Lcsr_df <- as.data.frame(L_childcare)
colour=c("#0D657D","#ee770d","#D3D3D3")
csr_plot <- ggplot(Lcsr_df, aes(r, obs-r))+
# plot observed value
geom_line(colour=c("#4d4d4d"))+
geom_line(aes(r,theo-r), colour="red", linetype = "dashed")+
# plot simulation envelopes
geom_ribbon(aes(ymin=lo-r,ymax=hi-r),alpha=0.1, colour=c("#91bfdb")) +
xlab("Distance r (m)") +
ylab("L(r)-r") +
geom_rug(data=Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,], sides="b", colour=colour[1]) +
geom_rug(data=Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,], sides="b", colour=colour[2]) +
geom_rug(data=Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,], sides="b", color=colour[3]) +
theme_tufte()+
ggtitle(title)
text1<-"Significant clustering"
text2<-"Significant segregation"
text3<-"Not significant clustering/segregation"
# the below conditional statement is required to ensure that the labels (text1/2/3) are assigned to the correct traces
if (nrow(Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,])==0){
if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text3, traces = 4) %>%
rangeslider()
}else if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text2, traces = 4) %>%
rangeslider()
}else {
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text2, traces = 4) %>%
style(text = text3, traces = 5) %>%
rangeslider()
}
} else if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text1, traces = 4) %>%
rangeslider()
} else{
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text1, traces = 4) %>%
style(text = text3, traces = 5) %>%
rangeslider()
}
} else{
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text1, traces = 4) %>%
style(text = text2, traces = 5) %>%
style(text = text3, traces = 6) %>%
rangeslider()
}